Apenas eco-regiões
%%HTML
<h1> Índice</h1>
<h3>
<a href="#1">clicar aqui!</a> para ir para as primeiras figuras
</h3><br>
<h3> outras figuras:</h3><br>
<a href="#varintermonbio">Variação inter-anual por mês em cada bioma</a>
<br>
<a href="#varintrabio"> Variação intra_anual em cada bioma</a>
# imports
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
import numpy as np
import matplotlib.cm as cm
import matplotlib.colors as cls
from pandas import Series, DataFrame
from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes, mark_inset
import matplotlib.ticker as ticker
from collections import defaultdict, namedtuple
from pandas import ExcelWriter
import xlsxwriter
import pandas as pd
import calendar
plt.rcParams["font.family"] = "Times New Roman"
def get_ecos():
m = Basemap(projection='cyl')
shpf = m.readshapefile('ecoregions2017/ecoregions2017','ecos',linewidth=0.1)
return m.ecos_info, m.ecos
months = "Jan Fev Mar Abr Mai Jun Jul Ago Set Out Nov Dez".split(' ')
ecos_info, ecos = get_ecos()
Eco = namedtuple('Eco', ['name', 'biome_name', 'biome_num', 'area'])
info_dict = {eco['ECO_ID']:Eco(eco['ECO_NAME'],eco['BIOME_NAME'],eco['BIOME_NUM'], eco['SHAPE_AREA']) for eco in ecos_info}
eco_names = []
for i in range(847):
name = info_dict[i].name
if type(name) == bytes:
name = name.decode('latin-1')
eco_names.append(name)
biome_names = []
for i in range(847):
biome_names.append(info_dict[i].biome_name)
biome_nums = []
for i in range(847):
biome_nums.append(info_dict[i].biome_num)
areas = []
for i in range(847):
areas.append(info_dict[i].area)
biome_nums[0] = 0.0
result = {}
for i in range(1,18):
suf = '20{0:02d}'.format(i)
in_file= 'result' + suf
result[suf] = DataFrame(np.load(in_file), columns =months)
for year in result:
result[year].index.name ='Eco_Id'
for year in result:
result[year] = result[year].transpose()
head = DataFrame()
head = head.assign(ECO_NAME = eco_names, BIOME_NUM=biome_nums,BIOME_NAME=biome_names,AREAS= areas).transpose()
result['Atributos'] = head
df_finall = pd.concat(result)
df_finall
df_finall.loc['Atributos'].loc['BIOME_NUM']
df_finall[:][:12*17]
biome_grouped = df_finall[:][:12*17].groupby(df_finall.loc['Atributos'].loc['BIOME_NUM'], axis=1)
biome_grouped.sum()
#writer = ExcelWriter('results_newstyle2.xlsx',engine='xlsxwriter')
#df_finall.to_excel(writer,sheet_name='Burned Pixels',engine='xlsxwriter')
#writer.save()
total_biome = biome_grouped.sum().sum()
total_biome
total_eco = df_finall.iloc[:12*17].sum()
total_eco
total_eco_area = total_eco/(areas)
df_finall[:][12*17:-1:2]
df_finall[:][12*17 +2: 12*17 +3]
df_finall.loc['Atributos'].loc['AREAS']
df_finall.loc['Atributos'].loc['AREAS'].groupby(biome_nums).sum()
df_finall.T[df_finall.T['Atributos']['BIOME_NUM'] == 14]['Atributos']['AREAS'].sum()
areas_biome = df_finall.loc['Atributos'].loc['AREAS'].groupby(biome_nums).sum()
areas_biome
total_biome_area = total_biome/areas_biome
total_biome.iloc[14]/areas_biome[14]
302628.0/27.665056
df_finall[1: 12*17:12]
months
monthly = [ df_finall[n: 12*17:12].sum()/areas for n in range(12)]
monthly
dict_monthly = dict(zip(months,monthly))
dict_monthly['Jan']
monthly_biome = [ df_finall[n: 12*17:12].sum() for n in range(12)]
monthly_biome = [mb.groupby(biome_nums).sum()/areas_biome for mb in monthly_biome]
monthly_biome
dict_monthly_biome = dict(zip(months,monthly_biome))
annual = [ df_finall[n*12: (n+1) *12].sum()/areas for n in range(17)]
annual
years = [str(y) for y in range(2001,2018)]
years
dict_annual = dict(zip(years,annual))
dict_annual['2002']
annual_biome = [ df_finall[n*12: (n+1) *12].sum() for n in range(17)]
annual_biome = [ ab.groupby(biome_nums).sum()/areas_biome for ab in annual_biome]
dict_annual_biome = dict(zip(years,annual_biome))
annual_total = [ df_finall[n*12: (n+1) *12].sum().sum() for n in range(17)]
fig, ax = plt.subplots()
ax.plot(years, annual_total)
labels = ax.get_xticklabels()
plt.tight_layout()
plt.setp(labels, rotation=30, fontsize=10);plt.show()
monthly_total = [ df_finall[n: 12*17:12].sum().sum() for n in range(12)]
fig, ax = plt.subplots()
#ax[0].plot( np.arange(1,13),monthly_total)
#labels = ax.get_xticklabels()
#plt.setp(labels, rotation=30, fontsize=10)
#plt.xticks(months,months)
ax.plot( np.arange(12),monthly_total)
ax.set_xticks(np.arange(12))
ax.set_xticklabels(months)
plt.show()
months
all_months = [df_finall.iloc[m].sum() for m in range(17*12)]
import itertools
ym = itertools.product(years,months)
ymlabels = [y + m for y, m in ym]
ymlabels
from statistics import mean
media = mean(all_months)
fig, ax = plt.subplots()
ax.plot( np.arange(12*17),all_months)
ax.axhline(y=media, color='r', linestyle='--')
ax.set_xticks(np.arange(12*17))
ax.set_xticklabels(ymlabels)
plt.setp(labels, rotation=60);
plt.show()
df_finall
annual_total
n = 0
df_finall[n*12: (n+1) *12].sum().sum()
ar_annual_biome = np.array(annual_biome)
#testing
lixo = np.array(annual_biome)
lixo.shape
ar_annual_biome.reshape(17,15)
anbio = ar_annual_biome.T
anbio[0]
annual_biome
fig, ax = plt.subplots(5,3, figsize =(18,12))
for i in range(5):
for j in range(3):
ax[i][j].plot(years, anbio[i*3 + j])
labels = ax[i][j].get_xticklabels()
plt.setp(labels, rotation=30, fontsize=10)
plt.show()
basket = 'BasketSemNomes.csv'
basket_df = pd.read_csv(basket,sep=';')
basket_df
basket_df.describe()
!cat ITDcluster.r
allpixels = np.load('../MODIS/MCD64A1/allpixels.npy')
allpixels[0]
annual_total[0]
an_total = np.array(annual_total[:-1])
diffs = (allpixels - an_total)*100/allpixels
diffs.max()
fig, ax = plt.subplots(3,1)
ax[0].plot(years[:-1], annual_total[:-1])
ax[0].set_title("Pixeis havendo selecção p/ ecoregião")
ax[1].plot(years[:-1], allpixels)
ax[1].set_title("Pixeis não havendo selecção p/ ecoregião")
ax[2].plot(years[:-1], diffs)
ax[2].set_title("Diferença em percentagem %")
labels = [ax[i].get_xticklabels() for i in range(3)]
plt.setp(labels, rotation=30, fontsize=10)
fig.tight_layout()
plt.show()
biome_names_dict = dict( (num, name) for num, name in zip(biome_nums, biome_names))
biome_names_dict
fig, ax = plt.subplots(5,3, figsize =(18,12))
for i in range(5):
for j in range(3):
ax[i][j].plot(years, anbio[i*3 + j])
labels = ax[i][j].get_xticklabels()
ax[i][j].set_title(biome_names_dict[i*3 + j])
plt.setp(labels, rotation=30, fontsize=10)
fig.tight_layout()
plt.show()
experimentar aplicar a escala das savanas às florestas mediterrânicas
savanas = 7; mediterraneo = 12
fig, ax = plt.subplots(2,1, sharey=True)
ax[0].plot(years, anbio[savanas])
labels = ax[0].get_xticklabels()
ax[0].set_title(biome_names_dict[savanas])
plt.setp(labels, rotation=30, fontsize=10)
ax[1].plot(years, anbio[mediterraneo])
labels = ax[1].get_xticklabels()
ax[1].set_title(biome_names_dict[mediterraneo])
plt.setp(labels, rotation=30, fontsize=10)
fig.tight_layout()
plt.show()
fig, ax = plt.subplots(2,1, sharey=True)
ver também sharey='col' e sharey='row'
log_anbio = (np.log(anbio) - np.log(anbio.min()))/(np.log(anbio.max()) - np.log(anbio.min()))
fig, ax = plt.subplots(2,1, sharey=True)
fig.suptitle('escala logarítmica normalizada')
ax[0].plot(years, log_anbio[savanas])
labels = ax[0].get_xticklabels()
ax[0].set_title(biome_names_dict[savanas])
plt.setp(labels, rotation=30, fontsize=10)
ax[1].plot(years, log_anbio[mediterraneo])
labels = ax[1].get_xticklabels()
ax[1].set_title(biome_names_dict[mediterraneo])
plt.setp(labels, rotation=30, fontsize=10)
fig.tight_layout()
plt.show()
log2_anbio = np.log(anbio -anbio.min() + 1) / np.log(anbio.max() - anbio.min() +1 )
fig, ax = plt.subplots(2,1, sharey=True)
fig.suptitle('escala logarítmica normalizada 2')
ax[0].plot(years, log2_anbio[savanas])
labels = ax[0].get_xticklabels()
ax[0].set_title(biome_names_dict[savanas])
plt.setp(labels, rotation=30, fontsize=10)
ax[1].plot(years, log2_anbio[mediterraneo])
labels = ax[1].get_xticklabels()
ax[1].set_title(biome_names_dict[mediterraneo])
plt.setp(labels, rotation=30, fontsize=10)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
x = np.linspace(.1,100000, 1000000)
log1 = np.log(x - x.min() +1) /np.log(x.max() - x.min() +1)
log2 = (np.log(x) - np.log(x.min()))/(np.log(x.max()) - np.log(x.min()))
fig, ax = plt.subplots(2,1, sharey=True)
fig.suptitle('escalas logarítmicas normalizadas')
ax[0].plot(x, log2, 'b')
ax[0].set_title(r"$\frac{\log(x) - \log(x_{min})} {\log(x_{max}) - \log(x_{min})}$" )
ax[1].plot(x, log1, 'r')
ax[1].set_title(r"$\frac{\log(x - x_{min} +1)}{ \log(x_{max}- x_{min} +1)}$")
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
fig, ax = plt.subplots(2,1, sharey=True,sharex=True )
labels = years[::2]
ax[0].plot(years, anbio[savanas])
ax[0].set_title(biome_names_dict[savanas])
ax[1].plot(years, anbio[mediterraneo])
ax[1].set_title(biome_names_dict[mediterraneo])
plt.xticks(labels,labels)
fig.tight_layout()
plt.show()
fig, ax = plt.subplots(2,1, sharey=True)
labels = years[::2]
#plt.setp(ax, xticks=labels, xticklabels=labels)
ax[0].plot(years, anbio[savanas])
ax[0].set_title(biome_names_dict[savanas])
ax[0].set_xticks(labels)
ax[0].set_xticklabels(labels)
ax[1].plot(years, anbio[mediterraneo])
ax[1].set_title(biome_names_dict[mediterraneo])
ax[1].set_xticks(labels)
ax[1].set_xticklabels(labels)
fig.tight_layout()
plt.show()
from scipy import stats
years_float = np.array([float(yr) for yr in years])
fig, ax = plt.subplots(5,3, figsize =(18,12) , sharey=True)
labels = years[::2]
for i in range(5):
for j in range(3):
slope, intercept, r_value, p_value, std_err = stats.linregress(years_float, anbio[i*3 + j])
ax[i][j].plot(years, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i][j].plot(years, anbio[i*3 + j], 'o')
ax[i][j].set_xticks(labels)
ax[i][j].set_xticklabels(labels)
ax[i][j].set_title(biome_names_dict[i*3 + j])
ax[i][j].legend()
fig.tight_layout()
plt.show()
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html
import seaborn as sns; sns.set()
fig, ax = plt.subplots(2,1, sharey=True)
labels = years[::2]
#plt.setp(ax, xticks=labels, xticklabels=labels)
ax[0].plot(years, anbio[savanas])
ax[0].set_title(biome_names_dict[savanas])
ax[0].set_xticks(labels)
ax[0].set_xticklabels(labels)
ax[1].plot(years, anbio[mediterraneo])
ax[1].set_title(biome_names_dict[mediterraneo])
ax[1].set_xticks(labels)
ax[1].set_xticklabels(labels)
fig.tight_layout()
plt.show()
fig, ax = plt.subplots(5,3, figsize =(18,12) , sharey=True)
labels = years[::2]
for i in range(5):
for j in range(3):
slope, intercept, r_value, p_value, std_err = stats.linregress(years_float, anbio[i*3 + j])
ax[i][j].plot(years, intercept + slope*years_float, 'r', label="r-squared: {} ".format(r_value**2))
ax[i][j].plot(years, anbio[i*3 + j], 'o')
ax[i][j].set_xticks(labels)
ax[i][j].set_xticklabels(labels)
ax[i][j].set_title(biome_names_dict[i*3 + j])
ax[i][j].legend()
fig.tight_layout()
plt.show()
%%HTML
<h1 id="1">PLOTS<h1>
fig, ax = plt.subplots()
labels = years[::2]
ax.plot(years, annual_total)
ax.set_xticks(labels)
ax.set_xticklabels(labels, fontsize=15)
ax.set_title('Interannual Variability',fontsize=15)
labels_x = ax.get_xticklabels()
labels_y = ax.get_yticklabels()
off =ax.get_yaxis().get_offset_text()
off.set_x(-.07)
off.set_fontsize(12)
#plt.setp(labels_x,fontsize=15, family='Times New Roman')
plt.setp(labels_y,fontsize=15)
#plt.xlabel(r"Burned Pixels", fontsize = 12)
plt.ylabel(r"Burned Pixels", fontsize = 15)
plt.tight_layout()
plt.show()
#help(off)
months = 'Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec'.split(' ')
fig, ax = plt.subplots()
ax.plot( np.arange(12),monthly_total)
ax.set_xticks(np.arange(12))
ax.set_xticklabels(months, fontsize=15)
ax.set_title('Intra-annual Variability',fontsize=15)
labels_x = ax.get_xticklabels()
labels_y = ax.get_yticklabels()
off =ax.get_yaxis().get_offset_text()
off.set_x(-.07)
off.set_fontsize(12)
plt.setp(labels_y,fontsize=15)
plt.ylabel(r"Burned Pixels", fontsize = 15)
plt.tight_layout()
plt.show()
residuals_monthly_total = (np.array(monthly_total) - mean(monthly_total))/mean(monthly_total)
fig, ax = plt.subplots()
ax.plot( np.arange(12),residuals_monthly_total)
ax.set_xticks(np.arange(12))
ax.set_xticklabels(months, fontsize=15)
ax.set_title('Intra-annual Variability',fontsize=15)
labels_x = ax.get_xticklabels()
labels_y = ax.get_yticklabels()
off =ax.get_yaxis().get_offset_text()
off.set_x(-.07)
off.set_fontsize(12)
plt.setp(labels_y,fontsize=15)
plt.ylabel(r"$\frac{Burned Pixels - \langle Burned Pixels \rangle}{\langle Burned Pixels \rangle}$", fontsize = 15)
plt.tight_layout()
plt.show()
130070.096/3600
$\frac{Burned Pixels - \langle Burned Pixels \rangle}{\langle Burned Pixels \rangle}$
ym = itertools.product(years,months)
ymlabels = [ m + ' ' + y[-2:] for y, m in ym]
fig, ax = plt.subplots(figsize=(16,12))
labels = ymlabels[::6]
ax.plot( np.arange(12*17),all_months)
ax.axhline(y=mean(all_months), color='r', linestyle='--', label='average')
ax.set_xticks(np.arange(0,12*17,6))
ax.set_xticklabels(labels,fontsize=12, rotation=30)
labels_y = ax.get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax.set_title('Intra-monthly Variability',fontsize=15)
plt.legend()
plt.show()
fig, ax = plt.subplots(5,3, figsize =(18,12) , sharey=True)
labels = years[::2]
for i in range(5):
for j in range(3):
slope, intercept, r_value, p_value, std_err = stats.linregress(years_float, anbio[i*3 + j])
ax[i][j].plot(years, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i][j].plot(years, anbio[i*3 + j], 'o')
ax[i][j].set_xticks(labels)
ax[i][j].set_xticklabels(labels,fontsize=15)
labels_y = ax[i][j].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i][j].set_title(biome_names_dict[i*3 + j],fontsize=15)
ax[i][j].legend(fontsize=15)
fig.tight_layout()
plt.show()
%%HTML
<h1 id="varintrabio"> Variação intra_anual em cada bioma</h1>
len(monthly_biome[0])
ar_monthly_biome = np.array(monthly_biome)
ar_monthly_biome.shape
monbio = ar_monthly_biome.T
monbio.shape
fig, ax = plt.subplots(5,3, figsize =(18,12) , sharey=True)
labels = months
for i in range(5):
for j in range(3):
#slope, intercept, r_value, p_value, std_err = stats.linregress(np.arange(12), monbio[i*3 + j])
#ax[i][j].plot(months, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i][j].plot(months, monbio[i*3 + j], 'o')
ax[i][j].set_xticks(np.arange(12))
ax[i][j].set_xticklabels(labels,fontsize=15)
labels_y = ax[i][j].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i][j].set_title(biome_names_dict[i*3 + j],fontsize=15)
ax[i][j].legend(fontsize=15)
fig.tight_layout()
plt.show()
fig, ax = plt.subplots(5,3, figsize =(18,12) )
labels = months
for i in range(5):
for j in range(3):
#slope, intercept, r_value, p_value, std_err = stats.linregress(np.arange(12), monbio[i*3 + j])
#ax[i][j].plot(months, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i][j].plot(months, monbio[i*3 + j], 'o')
ax[i][j].set_xticks(np.arange(12))
ax[i][j].set_xticklabels(labels,fontsize=15)
labels_y = ax[i][j].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i][j].set_title(biome_names_dict[i*3 + j],fontsize=15)
ax[i][j].legend(fontsize=15)
fig.tight_layout()
plt.show()
type(df_finall)
biome_grouped.sum()/areas_biome
biome_grouped.sum()[14][16*12]/areas_biome[14]
biomes_pa = biome_grouped.sum()/areas_biome
biomes_pa.index.get_level_values(level=1)
biomes_pa.index.get_level_values(level=1)
intra_biome = biomes_pa.groupby(biomes_pa.index.get_level_values(level=1), sort=False).sum()
for a, b in biomes_pa.groupby(biomes_pa.index.get_level_values(level=1),sort=False):
print(a)
print(b)
intra_biome
intra_biome[1]
intra_biome = [ b for a, b in biomes_pa.groupby(biomes_pa.index.get_level_values(level=1),sort=False)]
intra_biome
%%HTML
<h1 id="varintermonbio">Variação inter-anual por mês em cada bioma</h1>
labels = years[::2]
for b in range(15):
fig, ax = plt.subplots(4,3, figsize =(20,12) , sharey= True)
fig.suptitle(biome_names_dict[b],fontsize=15)
for i in range(4):
for j in range(3):
data = intra_biome[i*3 + j][b]
slope, intercept, r_value, p_value, std_err = stats.linregress(years_float, data)
ax[i][j].plot(years, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i][j].plot(years, data, 'o', label='_nolegend_')
ax[i][j].set_xticks(labels)
ax[i][j].set_xticklabels(labels,fontsize=15)
labels_y = ax[i][j].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i][j].set_title(calendar.month_name[i*3 + j +1],fontsize=15)
ax[i][j].legend(fontsize=15)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
calendar.month_name[3]
sns.set_context("paper")
labels = years[::2]
for b in range(15):
fig, ax = plt.subplots(4,3, figsize =(20,12) , sharey= True)
fig.suptitle(biome_names_dict[b],fontsize=15)
for i in range(4):
for j in range(3):
data = intra_biome[i*3 + j][b]
slope, intercept, r_value, p_value, std_err = stats.linregress(years_float, data)
ax[i][j].plot(years, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i][j].plot(years, data, 'o', label='_nolegend_')
ax[i][j].set_xticks(labels)
ax[i][j].set_xticklabels(labels,fontsize=15)
labels_y = ax[i][j].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i][j].set_title(calendar.month_name[i*3 + j +1],fontsize=15)
ax[i][j].legend(fontsize=15)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
sns.set_style("whitegrid")
for b in range(15):
fig, ax = plt.subplots(4,3, figsize =(20,12) , sharey= True)
fig.suptitle(biome_names_dict[b],fontsize=15)
for i in range(4):
for j in range(3):
data = intra_biome[i*3 + j][b]
slope, intercept, r_value, p_value, std_err = stats.linregress(years_float, data)
ax[i][j].plot(years, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
ax[i][j].plot(years, data, 'o', label='_nolegend_')
ax[i][j].set_xticks(labels)
ax[i][j].set_xticklabels(labels,fontsize=15)
labels_y = ax[i][j].get_yticklabels()
plt.setp(labels_y,fontsize=15)
ax[i][j].set_title(calendar.month_name[i*3 + j +1],fontsize=15)
ax[i][j].legend(fontsize=15)
fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
labels = years[::2]
for b in range(15):
fig, ax = plt.subplots(4,3, figsize =(20,12) , sharey= True)
fig.suptitle(biome_names_dict[b],fontsize=15)
for i in range(4):
for j in range(3):
data = intra_biome[i*3 + j][b]
#slope, intercept, r_value, p_value, std_err = stats.linregress(years_float, data)
#ax[i][j].plot(years, intercept + slope*years_float, 'r', label=r"$R^2$: {:5.3f} ".format(r_value**2))
sns.regplot(years_float, data, ax=ax[i][j])
#ax[i][j].set_xticks(labels)
#ax[i][j].set_xticklabels(labels,fontsize=15)
#labels_y = ax[i][j].get_yticklabels()
#plt.setp(labels_y,fontsize=15)
ax[i][j].set_title(calendar.month_name[i*3 + j +1],fontsize=15)
ax[i][j].set_ylabel('')
#ax[i][j].legend(fontsize=15)
#fig.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
x = np.linspace(1,200,10)
y = 4*x
df = pd.DataFrame(y, index=x)
df
df.columns = ['y']
df.index.name= 'x'
df
df =pd.DataFrame(dict(zip(['x','y'],[x,y])))
df
ax= sns.regplot( x = 'x', y='y',data=df);ax.set(ylabel=''); plt.show()
pwd
cd ../MODIS/MCD64A1
import datetime
from pyhdf.SD import SD, SDC
#import matplotlib.path as mplPath
import subprocess
import re
DATADIR = 'C:/Users/alpha/Documents/dados/MCD64A1/'
NECOS = 847
def get_tilehv(fn):
pattern = 'MCD64A1\.A\d{4}\d{3}\.h(\d{2})v(\d{2}).{22}'
# MCD64A1.A2005032.h14v04.006.2017017093958.hdf
matches= re.findall(pattern, fn)
if matches:
h, v = matches[0]
return int(h), int(v)
def get_coord(h,v):
lon_f = 'lon_coordh{}v{}.npy'.format(h,v)
lat_f = 'lat_coordh{}v{}.npy'.format(h,v)
lon = np.load(lon_f)
lat = np.load(lat_f)
return lon, lat
def get_julian(yr, mo):
'''
Return the first julian day of given month and year.
'''
dt = datetime.datetime(yr, mo, 1)
tt = dt.timetuple()
return tt.tm_yday
def julian2month(yr,jul):
# month from 1 to 12
values = [get_julian(yr,i) for i in range(1,13)]
#return month from 0 to 11
idx = [i for i in range(12) if values[i] == jul]
return idx[0]
def get_month(fn):
pattern = 'MCD64A1\.A(\d{4})(\d{3}).{29}'
matches= re.findall(pattern, fn)
year, julian = matches[0]
year, julian = int(year), int(julian)
return julian2month(year,julian)
file_ = DATADIR + 'MCD64A1.A2001152.h17v05.006.2017012103748.hdf'
h , v = get_tilehv(file_)
print('h: {}; v: {}'.format(h,v))
lon, lat = get_coord(h,v)
fich = SD(file_, SDC.READ)
# select sds
sds_obj = fich.select('Burn Date')
# get sds data
data = sds_obj.get()
attrs = sds_obj.attributes(full=1)
vra=attrs["valid_range"]
valid_range = vra[0]
fva=attrs["_FillValue"]
_FillValue = fva[0]
lat = lat.reshape(data.shape)
lon = lon.reshape(data.shape)
# Mask the data. logical_or from numpy
data = data.astype('double')
#invalid = logical_or(data < valid_range[0], data > valid_range[1])
invalid = np.logical_or(data < 1, data > valid_range[1])
invalid = np.logical_or(invalid, data == _FillValue)
# nan from numpy
data[invalid] = np.nan
# ma, isnan from numpy
data = np.ma.masked_array(data, np.isnan(data))
idx =np.logical_not(data.mask)
lat_val =lat[idx]
lon_val = lon[idx]
data_val = data[idx]
result = data_val.size
print('Result: {}'.format(result))
import pickle
from mpl_toolkits.basemap import Basemap
import matplotlib.path as mplPath
m = Basemap(projection='cyl')
shpf = m.readshapefile('ecoregions2017/ecoregions2017','ecos')
rings =[mplPath.Path(eco,closed=True) for eco in m.ecos]
ecos_info =m.ecos_info
result = np.zeros((NECOS),dtype=int)
br_vh = []
with open('bounded_rings.txt', 'rb') as f:
br_string = f.read()
br = pickle.loads(br_string)
br_vh = br[v,h]
not_included = []
for i, px in enumerate(data_val):
pt = ( lon_val[i], lat_val[i])
included = False
for ring_idx in br_vh:
ring = rings[ring_idx]
if ring.contains_point(pt):
ecoid = int(ecos_info[ring_idx]['ECO_ID'])
result[ecoid] += 1
included = True
break
if not included:
not_included.append(pt)
result.sum()
len(not_included)
total_ecoreg = df_finall.iloc[:12*17].sum() / df_finall.iloc[12*17]
total_ecoreg
ecos = m.ecos
cd -
fig, ax = plt.subplots(figsize=(12, 8), frameon= False)
for spine in plt.gca().spines.values():
spine.set_visible(False)
cmap = plt.get_cmap('OrRd')
min_val = 0
max_val = 1
norm = cls.Normalize(min_val, max_val)
cmmapable = cm.ScalarMappable(norm, cmap)
cmmapable.set_array(range(min_val, max_val))
ecos_info, ecos = get_ecos()
total_trans = np.log1p(total_ecoreg.values.astype('float'))/np.log1p(total_ecoreg.max())
patches = []
for info, ec in zip(ecos_info,ecos):
patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans[info['ECO_ID']] ))))
ax.add_collection(PatchCollection(patches, match_original= True, zorder=2))
title = 'Index Burned Area 2001-2017'
cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig('totalecologstd', dpi=1200, bbox_inches='tight')
#total_ecoreg.values.astype('float')
#np.log1p(total_ecoreg.values.astype('float'))
#total_ecoreg.values
pwd
np.log1p(total_ecoreg.max())
total_ecoreg.max()
total_ecoreg.size
areas_biome
df_finall.iloc[:12*17].groupby(df_finall.loc['Atributos'].loc['BIOME_NUM'], axis=1).sum().sum()/areas_biome
total_biome_area
fig, ax = plt.subplots(figsize=(12, 8), frameon= False)
for spine in plt.gca().spines.values():
spine.set_visible(False)
cmap = plt.get_cmap('OrRd')
min_val = 0
max_val = int(total_biome_area.max()) +1
norm = cls.Normalize(min_val, max_val)
cmmapable = cm.ScalarMappable(norm, cmap)
cmmapable.set_array(range(min_val, max_val))
ecos_info, ecos = get_ecos()
#total_trans = np.log(total_pa + 1)/np.log(total_pa.max() + 1)
patches = []
for info, ec in zip(ecos_info,ecos):
patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_biome_area[info['BIOME_NUM']] )))) if info['BIOME_NAME'] \
!= 'N/A' else patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_biome_area[0] ))))
ax.add_collection(PatchCollection(patches, match_original= True, zorder=2))
title = 'Biomes Burned Px / Area 2001-2017'
cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig('biomesnotransform2', dpi=1200, bbox_inches='tight')
df_finall[:][:12*17].T
df_finall.T['Atributos']['BIOME_NUM']
lala = pd.concat([df_finall[:][:12*17].T,df_finall.T['Atributos']['BIOME_NUM']], axis=1)
start = datetime.datetime.now()
fig, ax = plt.subplots(figsize=(12, 8), frameon= False)
for spine in plt.gca().spines.values():
spine.set_visible(False)
cmap = plt.get_cmap('OrRd')
min_val = 0
max_val = 1
norm = cls.Normalize(min_val, max_val)
cmmapable = cm.ScalarMappable(norm, cmap)
cmmapable.set_array(range(min_val, max_val))
ecos_info, ecos = get_ecos()
total_trans = np.log1p(total_biome_area)/np.log1p(total_biome_area.max())
patches = []
for info, ec in zip(ecos_info,ecos):
patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans[info['BIOME_NUM']] )))) if info['BIOME_NAME'] \
!= 'N/A' else patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans[0] )))) #cuidado estou a eliminar rock
ax.add_collection(PatchCollection(patches, match_original= True, zorder=2))
title = 'Biomes Index Burned 2001-2017'
cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig('biomeslogtransform2', dpi=1200, bbox_inches='tight')
print('It took {} minutes.'.format((datetime.datetime.now()-start)/60))
total_trans
import time
start = datetime.datetime.now()
time.sleep(29)
datetime.datetime.now()-start
print(str(datetime.timedelta(0, 29, 10693)).format('%m:%s'))
datetime.timedelta(0, 29, 10693).seconds
len([ df_finall[n*12: (n+1) *12].sum() for n in range(17)])
n= 1
df_finall[n*12: (n+1) *12].sum()
len(areas)
intereco =[ df_finall[n*12: (n+1) *12].sum()/areas for n in range(17)]
intereco_max= max(map(max,intereco))
for yr in range(17):
fig, ax = plt.subplots(figsize=(12, 8), frameon= False)
for spine in plt.gca().spines.values():
spine.set_visible(False)
cmap = plt.get_cmap('OrRd')
min_val = 0
max_val = 1
norm = cls.Normalize(min_val, max_val)
cmmapable = cm.ScalarMappable(norm, cmap)
cmmapable.set_array(range(min_val, max_val))
ecos_info, ecos = get_ecos()
total_trans = np.log1p(intereco[yr])/np.log1p(intereco_max)
patches = []
for info, ec in zip(ecos_info,ecos):
patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans[info['ECO_ID']] ))))
ax.add_collection(PatchCollection(patches, match_original= True, zorder=2))
title = 'Index Burned Area ' + years[yr]
cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig('ecologstd' + years[yr], dpi=1200, bbox_inches='tight')
annual_biome
type(annual_biome[0])
denominador = np.log1p(max(map(max,annual_biome)))
for yr in range(17):
fig, ax = plt.subplots(figsize=(12, 8), frameon= False)
for spine in plt.gca().spines.values():
spine.set_visible(False)
cmap = plt.get_cmap('OrRd')
min_val = 0
max_val = 1
norm = cls.Normalize(min_val, max_val)
cmmapable = cm.ScalarMappable(norm, cmap)
cmmapable.set_array(range(min_val, max_val))
ecos_info, ecos = get_ecos()
total_trans = np.log1p(annual_biome[yr])/denominador
patches = []
for info, ec in zip(ecos_info,ecos):
patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans.loc[info['BIOME_NUM']] )))) if info['BIOME_NAME'] \
!= 'N/A' else patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans.iloc[0] )))) #cuidado estou a eliminar rock
ax.add_collection(PatchCollection(patches, match_original= True, zorder=2))
title = 'Biomes Index Burned' + years[yr]
cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig('biomeslogtransform' + years[yr], dpi=1200, bbox_inches='tight')
annual_biome[0]
max(map(max,annual_biome))
open('df_finall','wb').write(pickle.dumps(df_finall))
lala = pickle.loads(open('df_finall','rb').read())
pwd
int(ecos_info[0]['BIOME_NUM'])
monthly_biome
notincluded2013 = pickle.loads(open('notincluded2013', 'rb').read())
notincluded2014 = pickle.loads(open('notincluded2014', 'rb').read())
len(notincluded2014)
lons, lats = zip(*notincluded2013)
lons=list(lons)
lats= list (lats)
plt.figure(figsize=(12,6))
map = Basemap(projection='cyl')
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()
x, y = map(lons, lats)
map.plot(lons,lats, 'mo', markersize=5)
plt.show()
lons, lats = zip(*notincluded2014)
plt.figure(figsize=(12,6))
map = Basemap(projection='cyl')
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()
x, y = map(lons, lats)
map.plot(lons,lats, 'mo', markersize=5)
plt.show()
lons14, lats14 = zip(*notincluded2014)
lons13, lats13 = zip(*notincluded2013)
plt.figure(figsize=(12,6))
map = Basemap(projection='cyl')
map.readshapefile('ecoregions2017/ecoregions2017','ecos',linewidth=0.1)
x, y = map(lons13, lats13)
map.plot(lons13,lats13, 'ro', markersize=2, label='2013')
map.plot(lons14,lats14, 'b+', markersize=2, label='2014')
plt.legend()
plt.show()
plt.figure(figsize=(12,6))
map = Basemap(projection='cyl')
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='white',lake_color='aqua')
map.drawcoastlines()
map.plot(lons13,lats13, 'go', markersize=2, label='2013')
map.plot(lons14,lats14, 'ro', markersize=2, label='2014')
plt.legend()
plt.show()
monthly_biome
max([ mb.max() for mb in monthly_biome])
np.log1p(20120.9686490036)
denominador2 = 9.9095674648050416
for mo in range(12):
fig, ax = plt.subplots(figsize=(12, 8), frameon= False)
for spine in plt.gca().spines.values():
spine.set_visible(False)
cmap = plt.get_cmap('OrRd')
min_val = 0
max_val = 1
norm = cls.Normalize(min_val, max_val)
cmmapable = cm.ScalarMappable(norm, cmap)
cmmapable.set_array(range(min_val, max_val))
ecos_info, ecos = get_ecos()
total_trans = np.log1p(monthly_biome[mo])/denominador2
patches = []
for info, ec in zip(ecos_info,ecos):
patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans.loc[info['BIOME_NUM']] )))) if info['BIOME_NAME'] \
!= 'N/A' else patches.append(Polygon(np.array(ec),True, color = cmap(norm(total_trans.iloc[0] )))) #cuidado estou a eliminar rock
ax.add_collection(PatchCollection(patches, match_original= True, zorder=2))
title = 'Biomes Index Burned ' + calendar.month_name[mo + 1]
cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig('biomeslogtransform' + calendar.month_name[mo + 1], dpi=1200, bbox_inches='tight')